Create a new analysis directories.
- general directory
- stains directory
- for plots
- for output of summary results
- for baseline tables
- for genetic analyses
- for Cox regression results
* General packages...
trying URL 'https://bioconductor.org/packages/3.21/bioc/bin/macosx/big-sur-x86_64/contrib/4.5/limma_3.64.1.tgz'
Content type 'application/x-gzip' length 3148172 bytes (3.0 MB)
==================================================
downloaded 3.0 MB
The downloaded binary packages are in
/var/folders/mm/b4fx9qss7t79pb95gn3k6grh0000gq/T//Rtmpp5gMOF/downloaded_packages
Update all/some/none? [a/s/n]:
For the ERA-CVD ‘druggable-MI-targets’ project (grantnumber: 01KL1802) we performed two related RNA sequencing (RNAseq) experiments:
conventional (‘bulk’) RNAseq using RNA extracted from carotid plaque samples, n ± 700. As of Wednesday, July 02, 2025 all samples have been selected and RNA has been extracted; quality control (QC) was performed and we have a dataset of 635 samples.
single-cell RNAseq (scRNAseq) of at least n = 40 samples (20 females, 20 males). As of Wednesday, July 02, 2025 data is available of 40 samples (3 females, 15 males), we are extending sampling to get more female samples.
Plaque samples are derived from carotid endarterectomies as part of the Athero-Express Biobank Study which is an ongoing study in the UMC Utrecht.
Here we map the PHENOMICL_downstream to single-cells from the plaques.
Here we obtain data from the PHENOMICL_downstream in plaques.
library(openxlsx)
gene_list_df <- read.xlsx(paste0(PROJECT_loc, "/targets/targets.xlsx"), sheet = "Genes")
gene_list <- unlist(gene_list_df$Gene)
gene_list [1] "SCGB3A2" "LIX1" "IGSF9B" "ND6" "ALB" "OR10A4" "RASEF" "NEDD4" "TRIM49D1"
[10] "TCL1A" "ND4L" "ATP8" "FBXO15" "F5" "TMEM212" "PTPRD" "CYP46A1" "OR2T33"
[19] "SORBS2" "ITGA7" "RNASE1" "FOS" "HMOX1" "LAPTM5" "MMP9" "SIGLEC1" "FTL"
[28] "CD14" "HCST" "TIMP3" "CCL2" "SAT1" "CD163" "PTGDS" "LGALS9" "ACKR1"
[37] "NT5DC2" "TGFBI" "C1QC" "S100A9"
First we will load the data:
Here we load the latest dataset from our Athero-Express single-cell RNA experiment.
# load(paste0(AESCRNA_loc, "/20210811.46.patients.KP.RData"))
# scRNAseqData <- seuset
# rm(seuset)
#
# saveRDS(scRNAseqData, paste0(AESCRNA_loc, "/20210811.46.patients.KP.RDS"))
scRNAseqData <- readRDS(paste0(AESCRNA_loc, "/20210811.46.patients.KP.RDS"))
scRNAseqDataAn object of class Seurat
36147 features across 4948 samples within 2 assays
Active assay: RNA (20111 features, 0 variable features)
3 layers present: counts, data, scale.data
1 other assay present: SCT
2 dimensional reductions calculated: pca, umap
The naming/classification is based on a combination conventional markers. We do not claim to know the exact identity of each cell, rather we refer to cells as ‘KIT+ Mast cells”-like cells. Likewise we refer to the cell clusters as ’communities’ of cells that exhibit similar properties, i.e. similar defining markers (e.g. KIT).
We will rename the cell types to human readable names.
### change names for clarity
backup.scRNAseqData = scRNAseqData
# get the old names to change to new names
UMAPPlot(scRNAseqData, label = FALSE, pt.size = 1.25, label.size = 4, group.by = "ident") [1] CD3+ T Cells I CD3+ T Cells IV
[3] CD34+ Endothelial Cells I CD3+ T Cells V
[5] CD3+CD56+ NK Cells II CD3+ T Cells VI
[7] CD68+IL18+TLR4+TREM2+ Resident macrophages CD3+CD56+ NK Cells I
[9] ACTA2+ Smooth Muscle Cells CD3+ T Cells II
[11] FOXP3+ T Cells CD34+ Endothelial Cells II
[13] CD3+ T Cells III CD68+CD1C+ Dendritic Cells
[15] CD68+CASP1+IL1B+SELL+ Inflammatory macrophages CD79A+ Class-switched Memory B Cells
[17] CD68+ABCA1+OLR1+TREM2+ Foam Cells CD68+KIT+ Mast Cells
[19] CD68+CD4+ Monocytes CD79+ Plasma B Cells
20 Levels: CD3+ T Cells I CD3+ T Cells II CD3+ T Cells III ... CD79+ Plasma B Cells
celltypes <- c("CD68+CD4+ Monocytes" = "CD68+CD4+ Mono",
"CD68+IL18+TLR4+TREM2+ Resident macrophages" = "CD68+IL18+TLR4+TREM2+ MRes",
"CD68+CD1C+ Dendritic Cells" = "CD68+CD1C+ DC",
"CD68+CASP1+IL1B+SELL+ Inflammatory macrophages" = "CD68+CASP1+IL1B+SELL MInf",
"CD68+ABCA1+OLR1+TREM2+ Foam Cells" = "CD68+ABCA1+OLR1+TREM2+ FC",
# T-cells
"CD3+ T Cells I" = "CD3+ TC I",
"CD3+ T Cells II" = "CD3+ TC II",
"CD3+ T Cells III" = "CD3+ TC III",
"CD3+ T Cells IV" = "CD3+ TC IV",
"CD3+ T Cells V" = "CD3+ TC V",
"CD3+ T Cells VI" = "CD3+ TC VI",
"FOXP3+ T Cells" = "FOXP3+ TC",
# Endothelial cells
"CD34+ Endothelial Cells I" = "CD34+ EC I",
"CD34+ Endothelial Cells II" = "CD34+ EC II",
# SMC
"ACTA2+ Smooth Muscle Cells" = "ACTA2+ SMC",
# NK Cells
"CD3+CD56+ NK Cells I" = "CD3+CD56+ NK I",
"CD3+CD56+ NK Cells II" = "CD3+CD56+ NK II",
# Mast
"CD68+KIT+ Mast Cells" = "CD68+KIT+ MC",
"CD79A+ Class-switched Memory B Cells" = "CD79A+ BCmem",
"CD79+ Plasma B Cells" = "CD79+ BCplasma")
scRNAseqData <- Seurat::RenameIdents(object = scRNAseqData,
celltypes)UMAPPlot(scRNAseqData, label = TRUE, pt.size = 1.25, label.size = 4, group.by = "ident",
repel = TRUE)Loading the Athero-Express clinical data.
# AEDB.CEA <- readRDS(file = paste0(OUT_loc, "/", Today,".",PROJECTNAME,".AEDB.CEA.RDS"))
AEDB.CEA <- readRDS(file = paste0(OUT_loc, "/20241219.",PROJECTNAME,".AEDB.CEA.RDS"))# Baseline table variables
basetable_vars = c("Hospital", "ORyear", "Artery_summary",
"Age", "Gender",
# "TC_finalCU", "LDL_finalCU", "HDL_finalCU", "TG_finalCU",
"TC_final", "LDL_final", "HDL_final", "TG_final",
# "hsCRP_plasma",
"systolic", "diastoli", "GFR_MDRD", "BMI",
"KDOQI", "BMI_WHO",
"SmokerStatus", "AlcoholUse",
"DiabetesStatus",
"Hypertension.selfreport", "Hypertension.selfreportdrug", "Hypertension.composite", "Hypertension.drugs",
"Med.anticoagulants", "Med.all.antiplatelet", "Med.Statin.LLD",
"Stroke_Dx", "sympt", "Symptoms.5G", "AsymptSympt", "AsymptSympt2G",
"Symptoms.Update2G", "Symptoms.Update3G", "indexsymptoms_latest_4g",
"restenos", "stenose",
"CAD_history", "PAOD", "Peripheral.interv",
"EP_composite", "EP_composite_time", "EP_major", "EP_major_time",
"MAC_rankNorm", "SMC_rankNorm", "Macrophages.bin", "SMC.bin",
"Neutrophils_rankNorm", "MastCells_rankNorm",
"IPH.bin", "VesselDensity_rankNorm",
"Calc.bin", "Collagen.bin",
"Fat.bin_10", "Fat.bin_40", "OverallPlaquePhenotype", "Plaque_Vulnerability_Index",
"PCSK9_plasma", "PCSK9_plasma_rankNorm")
basetable_bin = c("Gender", "Artery_summary",
"KDOQI", "BMI_WHO",
"SmokerStatus", "AlcoholUse",
"DiabetesStatus",
"Hypertension.selfreport", "Hypertension.selfreportdrug", "Hypertension.composite", "Hypertension.drugs",
"Med.anticoagulants", "Med.all.antiplatelet", "Med.Statin.LLD",
"Stroke_Dx", "sympt", "Symptoms.5G", "AsymptSympt", "AsymptSympt2G",
"Symptoms.Update2G", "Symptoms.Update3G", "indexsymptoms_latest_4g",
"restenos", "stenose",
"CAD_history", "PAOD", "Peripheral.interv",
"EP_major", "EP_composite", "Macrophages.bin", "SMC.bin",
"IPH.bin",
"Calc.bin", "Collagen.bin",
"Fat.bin_10", "Fat.bin_40", "OverallPlaquePhenotype", "Plaque_Vulnerability_Index")
# basetable_bin
basetable_con = basetable_vars[!basetable_vars %in% basetable_bin]
# basetable_conmetadata <- scRNAseqData@meta.data %>% as_tibble() %>% separate(orig.ident, c("Patient", NA))
scRNAseqDataMeta <- metadata %>% distinct(Patient, .keep_all = TRUE)
scRNAseqDataMetaAE <- merge(scRNAseqDataMeta, AEDB.CEA, by.x = "Patient", by.y = "STUDY_NUMBER", sort = FALSE, all.x = TRUE)
dim(scRNAseqDataMetaAE)[1] 46 1237
# Replace missing data
# Ref: https://cran.r-project.org/web/packages/naniar/vignettes/replace-with-na.html
require(naniar)
na_strings <- c("NA", "N A", "N / A", "N/A", "N/ A",
"Not Available", "Not available",
"missing",
"-999", "-99",
"No data available/missing", "No data available/Missing")
# Then you write ~.x %in% na_strings - which reads as “does this value occur in the list of NA strings”.
scRNAseqDataMetaAE %>%
replace_with_na_all(condition = ~.x %in% na_strings)cat("====================================================================================================")====================================================================================================
SELECTION THE SHIZZLE
- sanity checking PRIOR to selection
library(data.table)
require(labelled)
ae.gender <- to_factor(scRNAseqDataMetaAE$Gender)
ae.hospital <- to_factor(scRNAseqDataMetaAE$Hospital)
table(ae.gender, ae.hospital, dnn = c("Sex", "Hospital"), useNA = "ifany") Hospital
Sex St. Antonius, Nieuwegein UMC Utrecht <NA>
female 0 18 0
male 0 26 0
<NA> 0 0 2
ae.artery <- to_factor(scRNAseqDataMetaAE$Artery_summary)
table(ae.artery, ae.gender, dnn = c("Sex", "Artery"), useNA = "ifany") Artery
Sex female male
No artery known (yet), no surgery (patient ill, died, exited study), re-numbered to AAA 0 0
carotid (left & right) 18 25
femoral/iliac (left, right or both sides) 0 0
other carotid arteries (common, external) 0 1
carotid bypass and injury (left, right or both sides) 0 0
aneurysmata (carotid & femoral) 0 0
aorta 0 0
other arteries (renal, popliteal, vertebral) 0 0
femoral bypass, angioseal and injury (left, right or both sides) 0 0
<NA> 0 0
Artery
Sex <NA>
No artery known (yet), no surgery (patient ill, died, exited study), re-numbered to AAA 0
carotid (left & right) 0
femoral/iliac (left, right or both sides) 0
other carotid arteries (common, external) 0
carotid bypass and injury (left, right or both sides) 0
aneurysmata (carotid & femoral) 0
aorta 0
other arteries (renal, popliteal, vertebral) 0
femoral bypass, angioseal and injury (left, right or both sides) 0
<NA> 2
ae.gender
ae.ic female
missing 0
no, died 0
yes 9
yes, health treatment when possible 5
yes, no health treatment 2
yes, no health treatment, no commercial business 1
yes, no tissue, no commerical business 0
yes, no tissue, no questionnaires, no medical info, no commercial business 0
yes, no questionnaires, no health treatment, no commercial business 0
yes, no questionnaires, health treatment when possible 0
yes, no tissue, no questionnaires, no health treatment, no commerical business 0
yes, no health treatment, no medical info, no commercial business 0
yes, no tissue, no questionnaires, no health treatment, no medical info, no commercial business 0
yes, no questionnaires, no health treatment 0
yes, no tissue, no health treatment 0
yes, no tissue, no questionnaires 0
yes, no tissue, health treatment when possible 0
yes, no tissue 0
yes, no commerical business 1
yes, health treatment when possible, no commercial business 0
yes, no medical info, no commercial business 0
yes, no questionnaires 0
yes, no tissue, no questionnaires, no health treatment, no medical info 0
yes, no tissue, no questionnaires, no health treatment, no commercial business 0
yes, no medical info 0
yes, no questionnaires, no commercial business 0
yes, no questionnaires, no health treatment, no medical info 0
yes, no questionnaires, health treatment when possible, no commercial business 0
yes, no health treatment, no medical info 0
no, doesn't want to 0
no, unable to sign 0
no, no reaction 0
no, lost 0
no, too old 0
yes, no medical info, health treatment when possible 0
no (never asked for IC because there was no tissue) 0
yes, no medical info, no commercial business, health treatment when possible 0
no, endpoint 0
wil niets invullen, wel alles gebruiken 0
second informed concents: yes, no commercial business 0
nooit geincludeerd 0
yes, not outside EU 0
yes, no DNA 0
<NA> 0
ae.gender
ae.ic male
missing 0
no, died 0
yes 14
yes, health treatment when possible 7
yes, no health treatment 2
yes, no health treatment, no commercial business 2
yes, no tissue, no commerical business 0
yes, no tissue, no questionnaires, no medical info, no commercial business 0
yes, no questionnaires, no health treatment, no commercial business 0
yes, no questionnaires, health treatment when possible 0
yes, no tissue, no questionnaires, no health treatment, no commerical business 0
yes, no health treatment, no medical info, no commercial business 0
yes, no tissue, no questionnaires, no health treatment, no medical info, no commercial business 0
yes, no questionnaires, no health treatment 0
yes, no tissue, no health treatment 0
yes, no tissue, no questionnaires 0
yes, no tissue, health treatment when possible 0
yes, no tissue 0
yes, no commerical business 1
yes, health treatment when possible, no commercial business 0
yes, no medical info, no commercial business 0
yes, no questionnaires 0
yes, no tissue, no questionnaires, no health treatment, no medical info 0
yes, no tissue, no questionnaires, no health treatment, no commercial business 0
yes, no medical info 0
yes, no questionnaires, no commercial business 0
yes, no questionnaires, no health treatment, no medical info 0
yes, no questionnaires, health treatment when possible, no commercial business 0
yes, no health treatment, no medical info 0
no, doesn't want to 0
no, unable to sign 0
no, no reaction 0
no, lost 0
no, too old 0
yes, no medical info, health treatment when possible 0
no (never asked for IC because there was no tissue) 0
yes, no medical info, no commercial business, health treatment when possible 0
no, endpoint 0
wil niets invullen, wel alles gebruiken 0
second informed concents: yes, no commercial business 0
nooit geincludeerd 0
yes, not outside EU 0
yes, no DNA 0
<NA> 0
ae.gender
ae.ic <NA>
missing 0
no, died 0
yes 0
yes, health treatment when possible 0
yes, no health treatment 0
yes, no health treatment, no commercial business 0
yes, no tissue, no commerical business 0
yes, no tissue, no questionnaires, no medical info, no commercial business 0
yes, no questionnaires, no health treatment, no commercial business 0
yes, no questionnaires, health treatment when possible 0
yes, no tissue, no questionnaires, no health treatment, no commerical business 0
yes, no health treatment, no medical info, no commercial business 0
yes, no tissue, no questionnaires, no health treatment, no medical info, no commercial business 0
yes, no questionnaires, no health treatment 0
yes, no tissue, no health treatment 0
yes, no tissue, no questionnaires 0
yes, no tissue, health treatment when possible 0
yes, no tissue 0
yes, no commerical business 0
yes, health treatment when possible, no commercial business 0
yes, no medical info, no commercial business 0
yes, no questionnaires 0
yes, no tissue, no questionnaires, no health treatment, no medical info 0
yes, no tissue, no questionnaires, no health treatment, no commercial business 0
yes, no medical info 0
yes, no questionnaires, no commercial business 0
yes, no questionnaires, no health treatment, no medical info 0
yes, no questionnaires, health treatment when possible, no commercial business 0
yes, no health treatment, no medical info 0
no, doesn't want to 0
no, unable to sign 0
no, no reaction 0
no, lost 0
no, too old 0
yes, no medical info, health treatment when possible 0
no (never asked for IC because there was no tissue) 0
yes, no medical info, no commercial business, health treatment when possible 0
no, endpoint 0
wil niets invullen, wel alles gebruiken 0
second informed concents: yes, no commercial business 0
nooit geincludeerd 0
yes, not outside EU 0
yes, no DNA 0
<NA> 2
rm(ae.gender, ae.hospital, ae.artery, ae.ic)
scRNAseqDataMetaAE.all <- subset(scRNAseqDataMetaAE,
(Artery_summary == "carotid (left & right)" | Artery_summary == "other carotid arteries (common, external)" ) & # we only want carotids
informedconsent != "missing" & # we are really strict in selecting based on 'informed consent'!
informedconsent != "no, died" &
informedconsent != "yes, no tissue, no commerical business" &
informedconsent != "yes, no tissue, no questionnaires, no medical info, no commercial business" &
informedconsent != "yes, no tissue, no questionnaires, no health treatment, no commerical business" &
informedconsent != "yes, no tissue, no questionnaires, no health treatment, no medical info, no commercial business" &
informedconsent != "yes, no tissue, no health treatment" &
informedconsent != "yes, no tissue, no questionnaires" &
informedconsent != "yes, no tissue, health treatment when possible" &
informedconsent != "yes, no tissue" &
informedconsent != "yes, no tissue, no questionnaires, no health treatment, no medical info" &
informedconsent != "yes, no tissue, no questionnaires, no health treatment, no commercial business" &
informedconsent != "no, doesn't want to" &
informedconsent != "no, unable to sign" &
informedconsent != "no, no reaction" &
informedconsent != "no, lost" &
informedconsent != "no, too old" &
informedconsent != "yes, no medical info, health treatment when possible" &
informedconsent != "no (never asked for IC because there was no tissue)" &
informedconsent != "no, endpoint" &
informedconsent != "nooit geincludeerd" &
informedconsent != "yes, no health treatment, no commercial business" & # IMPORTANT: since we are sharing with a commercial party
informedconsent != "yes, no tissue, no commerical business" &
informedconsent != "yes, no tissue, no questionnaires, no medical info, no commercial business" &
informedconsent != "yes, no questionnaires, no health treatment, no commercial business" &
informedconsent != "yes, no tissue, no questionnaires, no health treatment, no commerical business" &
informedconsent != "yes, no health treatment, no medical info, no commercial business" &
informedconsent != "yes, no tissue, no questionnaires, no health treatment, no medical info, no commercial business" &
informedconsent != "yes, no commerical business" &
informedconsent != "yes, health treatment when possible, no commercial business" &
informedconsent != "yes, no medical info, no commercial business" &
informedconsent != "yes, no tissue, no questionnaires, no health treatment, no commercial business" &
informedconsent != "yes, no questionnaires, no commercial business" &
informedconsent != "yes, no questionnaires, health treatment when possible, no commercial business" &
informedconsent != "second informed concents: yes, no commercial business")
# scRNAseqDataMetaAE.all[1:10, 1:10]
dim(scRNAseqDataMetaAE.all)[1] 39 1237
Showing the baseline table for the scRNAseq data in 39 CEA patients with informed consent.
===========================================================================================
CREATE BASELINE TABLE
# Create baseline tables
# http://rstudio-pubs-static.s3.amazonaws.com/13321_da314633db924dc78986a850813a50d5.html
scRNAseqDataMetaAE.all.tableOne = print(CreateTableOne(vars = basetable_vars,
factorVars = basetable_bin,
strata = "Gender",
data = scRNAseqDataMetaAE.all, includeNA = TRUE),
nonnormal = c(),
quote = FALSE, showAllLevels = TRUE,
format = "p",
contDigits = 3) Stratified by Gender
level female
n 16
Hospital (%) St. Antonius, Nieuwegein 0.0
UMC Utrecht 100.0
ORyear (%) No data available/missing 0.0
2002 0.0
2003 0.0
2004 0.0
2005 0.0
2006 0.0
2007 0.0
2008 0.0
2009 0.0
2010 0.0
2011 0.0
2012 0.0
2013 0.0
2014 0.0
2015 0.0
2016 0.0
2017 0.0
2018 25.0
2019 43.8
2020 25.0
2021 6.2
2022 0.0
Artery_summary (%) carotid (left & right) 100.0
other carotid arteries (common, external) 0.0
Age (mean (SD)) 70.500 (8.017)
Gender (%) female 100.0
male 0.0
TC_final (mean (SD)) 5.180 (1.628)
LDL_final (mean (SD)) 3.444 (1.076)
HDL_final (mean (SD)) 1.200 (0.153)
TG_final (mean (SD)) 2.111 (0.818)
systolic (mean (SD)) 152.938 (30.732)
diastoli (mean (SD)) 81.375 (19.145)
GFR_MDRD (mean (SD)) 69.529 (22.867)
BMI (mean (SD)) 25.920 (4.284)
KDOQI (%) Normal kidney function 12.5
CKD 2 (Mild) 50.0
CKD 3 (Moderate) 31.2
<NA> 6.2
BMI_WHO (%) Underweight 0.0
Normal 43.8
Overweight 31.2
Obese 18.8
<NA> 6.2
SmokerStatus (%) Current smoker 31.2
Ex-smoker 50.0
Never smoked 12.5
<NA> 6.2
AlcoholUse (%) No 43.8
Yes 50.0
<NA> 6.2
DiabetesStatus (%) Control (no Diabetes Dx/Med) 87.5
Diabetes 12.5
Hypertension.selfreport (%) no 0.0
yes 87.5
<NA> 12.5
Hypertension.selfreportdrug (%) no 0.0
yes 93.8
<NA> 6.2
Hypertension.composite (%) no 6.2
yes 93.8
Hypertension.drugs (%) no 12.5
yes 81.2
<NA> 6.2
Med.anticoagulants (%) no 87.5
yes 0.0
<NA> 12.5
Med.all.antiplatelet (%) no 18.8
yes 75.0
<NA> 6.2
Med.Statin.LLD (%) no 25.0
yes 68.8
<NA> 6.2
Stroke_Dx (%) No stroke diagnosed 56.2
Stroke diagnosed 43.8
sympt (%) Asymptomatic 0.0
TIA 31.2
minor stroke 25.0
Major stroke 12.5
Amaurosis fugax 12.5
Retinal infarction 6.2
Symptomatic, but aspecific symtoms 6.2
retinal infarction 0.0
Ocular ischemic syndrome 6.2
Symptoms.5G (%) Asymptomatic 0.0
Ocular 18.8
Other 6.2
Retinal infarction 6.2
Stroke 37.5
TIA 31.2
AsymptSympt (%) Asymptomatic 0.0
Ocular and others 31.2
Symptomatic 68.8
AsymptSympt2G (%) Asymptomatic 0.0
Symptomatic 100.0
Symptoms.Update2G (%) Asymptomatic 0.0
Symptomatic 100.0
Symptoms.Update3G (%) Asymptomatic 0.0
Symptomatic 100.0
indexsymptoms_latest_4g (%) asymp 0.0
ocular 25.0
TIA 37.5
stroke 37.5
restenos (%) de novo 100.0
stenose (%) 0-49% 6.2
50-70% 0.0
70-90% 62.5
90-99% 25.0
70-99% 6.2
CAD_history (%) No history CAD 81.2
History CAD 18.8
PAOD (%) no 81.2
yes 18.8
Peripheral.interv (%) no 75.0
yes 25.0
EP_composite (%) No composite endpoints 81.2
Composite endpoints 12.5
<NA> 6.2
EP_composite_time (mean (SD)) 1302.998 (4867.095)
EP_major (%) No major events (endpoints) 87.5
Major events (endpoints) 6.2
<NA> 6.2
EP_major_time (mean (SD)) 1303.167 (4867.046)
Macrophages.bin (%) no/minor 0.0
<NA> 100.0
SMC.bin (%) moderate/heavy 0.0
<NA> 100.0
Calc.bin (%) no/minor 0.0
<NA> 100.0
Collagen.bin (%) moderate/heavy 0.0
<NA> 100.0
Fat.bin_10 (%) >10% 0.0
<NA> 100.0
Fat.bin_40 (%) <40% 0.0
<NA> 100.0
OverallPlaquePhenotype (%) fibroatheromatous 0.0
<NA> 100.0
Plaque_Vulnerability_Index (%) 0 100.0
1 0.0
Stratified by Gender
male p test
n 23
Hospital (%) 0.0 NaN
100.0
ORyear (%) 0.0 NaN
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
69.6
30.4
0.0
0.0
0.0
Artery_summary (%) 95.7 1.000
4.3
Age (mean (SD)) 73.174 (8.294) 0.322
Gender (%) 0.0 <0.001
100.0
TC_final (mean (SD)) 4.153 (0.800) 0.037
LDL_final (mean (SD)) 2.244 (0.687) 0.002
HDL_final (mean (SD)) 1.094 (0.262) 0.257
TG_final (mean (SD)) 1.816 (1.244) 0.534
systolic (mean (SD)) 149.318 (24.602) 0.689
diastoli (mean (SD)) 78.500 (14.500) 0.601
GFR_MDRD (mean (SD)) 86.166 (36.012) 0.127
BMI (mean (SD)) 26.627 (3.794) 0.604
KDOQI (%) 39.1 0.161
21.7
26.1
13.0
BMI_WHO (%) 4.3 0.730
26.1
43.5
17.4
8.7
SmokerStatus (%) 26.1 0.972
56.5
13.0
4.3
AlcoholUse (%) 34.8 0.841
56.5
8.7
DiabetesStatus (%) 60.9 0.145
39.1
Hypertension.selfreport (%) 13.0 0.084
87.0
0.0
Hypertension.selfreportdrug (%) 13.0 0.319
82.6
4.3
Hypertension.composite (%) 8.7 1.000
91.3
Hypertension.drugs (%) 8.7 0.889
87.0
4.3
Med.anticoagulants (%) 87.0 0.332
8.7
4.3
Med.all.antiplatelet (%) 21.7 0.947
73.9
4.3
Med.Statin.LLD (%) 17.4 0.799
78.3
4.3
Stroke_Dx (%) 56.5 1.000
43.5
sympt (%) 26.1 0.215
8.7
26.1
8.7
17.4
0.0
0.0
4.3
8.7
Symptoms.5G (%) 26.1 0.126
26.1
0.0
4.3
34.8
8.7
AsymptSympt (%) 26.1 0.071
30.4
43.5
AsymptSympt2G (%) 26.1 0.077
73.9
Symptoms.Update2G (%) 26.1 0.077
73.9
Symptoms.Update3G (%) 26.1 0.077
73.9
indexsymptoms_latest_4g (%) 26.1 0.044
30.4
8.7
34.8
restenos (%) 100.0 NA
stenose (%) 0.0 0.122
17.4
34.8
26.1
21.7
CAD_history (%) 78.3 1.000
21.7
PAOD (%) 87.0 0.972
13.0
Peripheral.interv (%) 78.3 1.000
21.7
EP_composite (%) 82.6 0.965
13.0
4.3
EP_composite_time (mean (SD)) 2.330 (1.242) 0.226
EP_major (%) 87.0 0.932
8.7
4.3
EP_major_time (mean (SD)) 2.357 (1.254) 0.226
Macrophages.bin (%) 4.3 1.000
95.7
SMC.bin (%) 4.3 1.000
95.7
Calc.bin (%) 4.3 1.000
95.7
Collagen.bin (%) 4.3 1.000
95.7
Fat.bin_10 (%) 4.3 1.000
95.7
Fat.bin_40 (%) 4.3 1.000
95.7
OverallPlaquePhenotype (%) 4.3 1.000
95.7
Plaque_Vulnerability_Index (%) 95.7 1.000
4.3
Let’s save the baseline characteristics of the Athero-Express Biobank Study.
library(gtsummary)
library(dplyr)
library(haven)
# fix factor levels
scRNAseqDataMetaAE.all <- scRNAseqDataMetaAE.all %>%
mutate(across(where(haven::is.labelled), haven::as_factor))# Remove continuous variables with only one unique value (including NAs)
valid_vars <- scRNAseqDataMetaAE.all %>%
gtsummary::select(all_of(basetable_vars)) %>%
gtsummary::select(where(~n_distinct(., na.rm = TRUE) > 1)) %>%
names()
# Rebuild the table with valid variables
# These warnings:
# `cannot compute exact p-value with ties`
# are non-fatal and expected when using non-parametric tests (e.g., Wilcoxon) on ranked or integer-valued data that may contain duplicates (ties)
# To control this manually:
# add_p(
# test = all_continuous() ~ "t.test", # or "wilcox.test"
# test.args = all_categorical() ~ list(simulate.p.value = TRUE)
# )
# To have `gtsummary auto-select:
# add_p()
scRNAseqDataMetaAE.all.tbl <- scRNAseqDataMetaAE.all %>%
gtsummary::select(all_of(valid_vars), Gender) %>%
tbl_summary(
by = Gender,
missing = "ifany",
statistic = all_continuous() ~ "{mean} ({sd})"
) %>%
add_p(
test = all_continuous() ~ "t.test", # or "wilcox.test"
test.args = all_categorical() ~ list(simulate.p.value = TRUE)
) %>%
bold_labels()library(writexl)
# Convert to a clean data frame
scRNAseqDataMetaAE.all.tbl_df <- scRNAseqDataMetaAE.all.tbl %>% as_tibble()
# Save to Excel
write.xlsx(scRNAseqDataMetaAE.all.tbl_df,
file = paste0(BASELINE_loc, "/",Today,".",PROJECTNAME,".AESCRNA.CEA.39pts.after_qc.IC_academic.BaselineTable.xlsx"),
rowNames = FALSE,
colNames = TRUE,
sheetName = "AESCRNA_byGender", overwrite = TRUE)library(officer)
library(flextable)
# Format for Word output with Calibri font and size 9
ft <- gtsummary::as_flex_table(scRNAseqDataMetaAE.all.tbl) %>%
theme_vanilla() %>%
autofit() %>%
set_table_properties(layout = "autofit") %>%
flextable::font(fontname = "Calibri", part = "all") %>%
flextable::fontsize(size = 9, part = "all")
# Save to Word
doc <- read_docx() %>%
body_add_par("Baseline characteristics Athero-Express SCRNA Study -- stratified by Gender") %>%
body_add_flextable(ft)
print(doc, target = paste0(BASELINE_loc, "/",Today,".",PROJECTNAME,".AESCRNA.CEA.39pts.after_qc.IC_academic.BaselineTable.docx"))
# Print message
message("Table exported successfully to *.AESCRNA.CEA.39pts.after_qc.IC_academic.BaselineTable.xlsx and *.AESCRNA.CEA.39pts.after_qc.IC_academic.BaselineTable.docx with Calibri font and size 9.")
rm(ft, doc)List of samples to be included based on informed consent (see above).
[1] "4440" "4447" "4443" "4450" "4453" "4470" "4459" "4458" "4455" "4487" "4472" "4480" "4477" "4448"
[15] "4486" "4488" "4502" "4500" "4501" "4530" "4541" "4542" "4513" "4535" "4546" "4545" "4571" "4489"
[29] "4491" "4495" "4496" "4521" "4520" "4580" "4602" "4605" "4601" "4676" "4653"
Just make sure to keep the rows
(patientid.plateid.cellid) in the
data, so we will create a new column for that.
Subsetting the samples of interest based on
samples_of_interest.
SCTransform()Since we subsetted we need to re-run SCTransform().
# Extract the SCTModel
sct_model <- scRNAseqData@assays$SCT@SCTModel.list
# Inspect the regression variables
sct_model$model1
An sctransform model.
Model formula: y ~ log_umi
Parameters stored for 16036 features, 4948 cells.
[1] "RNA"
|
| | 0%
|
|======================= | 25%
|
|============================================== | 50%
|
|====================================================================== | 75%
|
|=============================================================================================| 100%
Selecting the clinical data of interest.
variables_of_interest <- c("Hospital", "ORyear", "Artery_summary",
"Age", "Gender",
"TC_final", "LDL_final", "HDL_final", "TG_final",
"systolic", "diastoli", "GFR_MDRD", "BMI",
"KDOQI", "BMI_WHO",
"SmokerStatus", "AlcoholUse",
"DiabetesStatus",
"Hypertension.selfreport", "Hypertension.selfreportdrug", "Hypertension.composite", "Hypertension.drugs",
"Med.anticoagulants", "Med.all.antiplatelet", "Med.Statin.LLD",
"Stroke_Dx",
"sympt", "Symptoms.5G", "AsymptSympt", "AsymptSympt2G",
"Symptoms.Update2G", "Symptoms.Update3G", "indexsymptoms_latest_4g",
"restenos", "stenose",
"CAD_history", "PAOD", "Peripheral.interv",
"EP_composite", "EP_composite_time", "EP_major", "EP_major_time")
temp <- subset(scRNAseqDataMetaAE.all, select = c("Patient", variables_of_interest))Now we are merging in the clinical data of interest.
# Step 1: Extract original metadata from Seurat object
original_metadata_df <- scRNAseqDataCEA39@meta.data
# Step 2: Merge with dataframe (temp) with select variables based on 'Patient' column
merged_metadata_df <- merge(original_metadata_df, temp, by = "Patient", all.x = TRUE, sort = FALSE)
# Step 3: Exclude already existing columns to avoid duplication
new_metadata_columns <- setdiff(colnames(temp), colnames(original_metadata_df))
new_meta_info_df <- merged_metadata_df[, new_metadata_columns, drop = FALSE]
# Step 4: Add the new metadata to the Seurat object
scRNAseqDataCEA39 <- AddMetaData(scRNAseqDataCEA39, metadata = new_meta_info_df)
# Step 5: Rename `Patient` to `STUDY_NUMBER` to match the clinical data
scRNAseqDataCEA39@meta.data <- dplyr::rename(scRNAseqDataCEA39@meta.data, "STUDY_NUMBER" = "Patient")
# Verify that the new metadata was added correctly
head(scRNAseqDataCEA39@meta.data)temp2 <- as_tibble(subset(scRNAseqDataCEA39@meta.data, select = c("STUDY_NUMBER", "orig.ident", "nCount_RNA", "nFeature_RNA",
"Plate", "Batch", "C.H", "Type", "percent.mt",
"nCount_SCT", "nFeature_SCT", "seurat_clusters")))
# fwrite(temp2,
# file = paste0(OUT_loc, "/", Today, ".AESCRNA.CEA.39pts.samplelist.after_qc.IC_commercial.csv"),
# sep = ",", row.names = FALSE, col.names = TRUE,
# showProgress = TRUE)
# rm(temp2)
#
# temp <- dplyr::rename(temp, "STUDY_NUMBER" = "Patient")
# fwrite(temp,
# file = paste0(OUT_loc, "/", Today, ".AESCRNA.CEA.39pts.clinicaldata.after_qc.IC_commercial.csv"),
# sep = ",", row.names = FALSE, col.names = TRUE,
# showProgress = TRUE)
# rm(temp)
#
# saveRDS(scRNAseqDataCEA39, file = paste0(OUT_loc, "/", Today, ".AESCRNA.CEA.39pts.Seurat.after_qc.IC_commercial.RDS"))
fwrite(temp2,
file = paste0(OUT_loc, "/", Today, ".AESCRNA.CEA.39pts.samplelist.after_qc.IC_academic.csv"),
sep = ",", row.names = FALSE, col.names = TRUE,
showProgress = TRUE)
rm(temp2)
temp <- dplyr::rename(temp, "STUDY_NUMBER" = "Patient")
fwrite(temp,
file = paste0(OUT_loc, "/", Today, ".AESCRNA.CEA.39pts.clinicaldata.after_qc.IC_academic.csv"),
sep = ",", row.names = FALSE, col.names = TRUE,
showProgress = TRUE)
rm(temp)
saveRDS(scRNAseqDataCEA39, file = paste0(OUT_loc, "/", Today, ".AESCRNA.CEA.39pts.Seurat.after_qc.IC_academic.RDS"))Here review the number of cells per sample, plate, and patients. And plot the ratio’s per sample and study number.
How many cells per type ...?
16 22 23 21 12 19 18 17 15 2 10 11 9 13 5 8 4 6 3 0 1
8 17 19 23 26 55 78 94 112 172 173 176 198 199 242 243 293 297 348 831 846
# cat("\n\nHow many cells per plate ...?")
# sort(table(scRNAseqDataCEA39@meta.data$ID))
# cat("\n\nHow many cells per type per plate ...?")
# table(scRNAseqDataCEA39@meta.data$SCT_snn_res.0.8, scRNAseqDataCEA39@meta.data$ID)
cat("\n\nHow many cells per patient ...?")
How many cells per patient ...?
4530 4440 4605 4653 4472 4458 4455 4496 4601 4502 4501 4571 4448 4477 4459 4520 4602 4489 4495 4545
3 6 7 20 22 35 54 70 70 73 75 76 80 84 94 96 96 97 102 106
4480 4447 4500 4513 4535 4676 4486 4470 4487 4546 4488 4521 4580 4491 4541 4450 4542 4453 4443
112 114 116 123 130 135 137 144 144 144 146 161 163 175 178 205 213 222 422
Visualizing these ratio's per study number and sample ...?
UMAPPlot(scRNAseqDataCEA39, label = TRUE, pt.size = 1.25, label.size = 4, group.by = "ident",
repel = TRUE)
ggsave(paste0(PLOT_loc, "/", Today, ".UMAP.png"), plot = last_plot())
# barplot(prop.table(x = table(scRNAseqDataCEA39@active.ident, scRNAseqDataCEA39@meta.data$Patient)),
# cex.axis = 1.0, cex.names = 0.5, las = 1,
# col = uithof_color, xlab = "study number", legend.text = FALSE, args.legend = list(x = "bottom"))
# dev.copy(pdf, paste0(QC_loc, "/", Today, ".cell_ratios_per_sample.pdf"))
# dev.off()
# barplot(prop.table(x = table(scRNAseqDataCEA39@active.ident, scRNAseqDataCEA39@meta.data$ID)),
# cex.axis = 1.0, cex.names = 0.5, las = 2,
# col = uithof_color, xlab = "sample ID", legend.text = FALSE, args.legend = list(x = "bottom"))
# dev.copy(pdf, paste0(QC_loc, "/", Today, ".cell_ratios_per_sample_per_plate.pdf"))
# dev.off()Let’s project known cellular markers.
UMAPPlot(scRNAseqDataCEA39, label = FALSE, pt.size = 1.25, label.size = 4, group.by = "ident",
repel = TRUE)
# endothelial cells
FeaturePlot(scRNAseqDataCEA39, features = c("CD34"), cols = c("#ECECEC", "#DB003F"))
# macrophages
FeaturePlot(scRNAseqDataCEA39, features = c("CD14", "CD68"), cols = c("#ECECEC", "#DB003F"))# FeaturePlot(scRNAseqDataCEA39, features = c("CD8"), cols = c("#ECECEC", "#DB003F"))
# b-cells
FeaturePlot(scRNAseqDataCEA39, features = c("CD79A"), cols = c("#ECECEC", "#DB003F"))We check whether the targets genes were sequenced using our method. In case some genes are not available in our data we could filter them here.
[1] "SCGB3A2" "LIX1" "IGSF9B" "ND6" "ALB" "OR10A4" "RASEF" "NEDD4" "TRIM49D1"
[10] "TCL1A" "ND4L" "ATP8" "FBXO15" "F5" "TMEM212" "PTPRD" "CYP46A1" "OR2T33"
[19] "SORBS2" "ITGA7" "RNASE1" "FOS" "HMOX1" "LAPTM5" "MMP9" "SIGLEC1" "FTL"
[28] "CD14" "HCST" "TIMP3" "CCL2" "SAT1" "CD163" "PTGDS" "LGALS9" "ACKR1"
[37] "NT5DC2" "TGFBI" "C1QC" "S100A9"
This code is just an example to filter the list from genes that are not in the data.
gene_list_rm <- c("ND6", "TRIM49D1",
"ND4L", "ATP8", "RNASE1")
temp = target_genes[!target_genes %in% gene_list_rm]
target_genes_qc <- c(temp)
# gene_list_qc <- gene_list
#
# for debug
# gene_list_qc_replace <- c("MRTFA")
# target_genes_qc <- target_genes
target_genes_qc [1] "SCGB3A2" "LIX1" "IGSF9B" "ALB" "OR10A4" "RASEF" "NEDD4" "TCL1A" "FBXO15"
[10] "F5" "TMEM212" "PTPRD" "CYP46A1" "OR2T33" "SORBS2" "ITGA7" "FOS" "HMOX1"
[19] "LAPTM5" "MMP9" "SIGLEC1" "FTL" "CD14" "HCST" "TIMP3" "CCL2" "SAT1"
[28] "CD163" "PTGDS" "LGALS9" "ACKR1" "NT5DC2" "TGFBI" "C1QC" "S100A9"
library(RColorBrewer)
p1 <- DotPlot(scRNAseqDataCEA39,
features = target_genes_qc,
cols = "RdBu") +
theme(axis.text.x = element_text(angle = 45, hjust=1, size = 8)) +
theme(axis.title.x = element_blank(), axis.title.y = element_blank()) # remove the y- and x-axis labels
p1
ggsave(paste0(PLOT_loc, "/", Today, ".DotPlot.Targets.png"), plot = last_plot())ggsave(paste0(PLOT_loc, "/", Today, ".DotPlot.Targets.ps"), plot = last_plot())
ggsave(paste0(PLOT_loc, "/", Today, ".DotPlot.Targets.pdf"), plot = last_plot())Let’s group this dotplot based on some grouping.
If we want we can also filter the genes that are not found in the data.
# Assuming 'gene_list_df' is your data.frame with the relevant columns loaded
gene_list_df_filtered <- subset(gene_list_df, Comments != "not found" | is.na(Comments))
gene_list_df_filtered# Separate genes by significance
gene_list_df_filtered_up <- gene_list_df_filtered$GeneSymbol[gene_list_df_filtered$significance == "UP"]
gene_list_df_filtered_down <- gene_list_df_filtered$GeneSymbol[gene_list_df_filtered$significance == "DOWN"]
# Combine the genes, first DOWN then UP (or vice versa based on preference)
ordered_genes_qc <- c(gene_list_df_filtered_up, gene_list_df_filtered_down)library(RColorBrewer)
# Create DotPlot ordered
p1_group <- DotPlot(scRNAseqData,
features = ordered_genes_qc,
cols = "RdBu") +
theme(axis.text.x = element_text(angle = 45, hjust=1, size = 8)) +
theme(axis.title.x = element_blank(), axis.title.y = element_blank()) # remove the y- and x-axis labels
p1_group
# Save the plots
ggsave(paste0(PLOT_loc, "/", Today, ".DotPlot.Targets.Grouped.png"), plot = p1_group)ggsave(paste0(PLOT_loc, "/", Today, ".DotPlot.Targets.Grouped.ps"), plot = p1_group)
ggsave(paste0(PLOT_loc, "/", Today, ".DotPlot.Targets.Grouped.pdf"), plot = p1_group)We should also plot these data stratified by gender smooth muscle cells.
[1] male female
Levels: female male
# Subset Seurat object for Males
male_cells <- WhichCells(scRNAseqDataCEA39, expression = Gender == "male")
scRNAseqDataCEA39_males <- subset(scRNAseqDataCEA39, cells = male_cells)
# Subset Seurat object for Females
female_cells <- WhichCells(scRNAseqDataCEA39, expression = Gender == "female")
scRNAseqDataCEA39_females <- subset(scRNAseqDataCEA39, cells = female_cells)Get an overview of the various cell communities.
[1] CD3+ TC I CD3+ TC IV CD34+ EC I
[4] CD3+ TC V CD3+CD56+ NK II CD3+ TC VI
[7] CD68+IL18+TLR4+TREM2+ MRes CD3+CD56+ NK I ACTA2+ SMC
[10] CD3+ TC II FOXP3+ TC CD34+ EC II
[13] CD3+ TC III CD68+CD1C+ DC CD68+CASP1+IL1B+SELL MInf
[16] CD79A+ BCmem CD68+ABCA1+OLR1+TREM2+ FC CD68+KIT+ MC
[19] CD68+CD4+ Mono CD79+ BCplasma
20 Levels: CD68+CD4+ Mono CD68+IL18+TLR4+TREM2+ MRes CD68+CD1C+ DC ... CD79+ BCplasma
Select the cell community of interest - in this example ‘smooth
muscle cells’ defined by ACTA2+ SMC.
# Define the cell identities you want to include in the plot
selected_idents_cells <- c("ACTA2+ SMC")library(RColorBrewer)
library(cowplot) # For combining plots
# Create DotPlot for males
p1_males_cells <- DotPlot(scRNAseqDataCEA39_males,
features = ordered_genes_qc,
idents = selected_idents_cells,
cols = "RdBu") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8)) +
theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
ggtitle("Males")
# Create DotPlot for females
p2_females_cells <- DotPlot(scRNAseqDataCEA39_females,
features = ordered_genes_qc,
idents = selected_idents_cells,
cols = "RdBu") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8)) +
theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
# theme(axis.title.x = element_blank(),
# axis.title.y = element_blank(),
# axis.text.y = element_blank(), # Remove y-axis text
# axis.ticks.y = element_blank()) + # Remove y-axis ticks
ggtitle("Females")
# Combine the two plots into one panel using cowplot
combined_plot_cells <- plot_grid(p1_males_cells, p2_females_cells, ncol = 2)
# Display the combined plot
print(combined_plot_cells)
# Save the plots
ggsave(paste0(PLOT_loc, "/", Today, ".DotPlot.Targets.ColocGrouped.TCellonly.GenderMales.png"), plot = p1_males_cells)ggsave(paste0(PLOT_loc, "/", Today, ".DotPlot.Targets.ColocGrouped.TCellonly.GenderMales.ps"), plot = p1_males_cells)
ggsave(paste0(PLOT_loc, "/", Today, ".DotPlot.Targets.ColocGrouped.TCellonly.GenderMales.pdf"), plot = p1_males_cells)
ggsave(paste0(PLOT_loc, "/", Today, ".DotPlot.Targets.ColocGrouped.TCellonly.GenderFemales.png"), plot = p2_females_cells)
ggsave(paste0(PLOT_loc, "/", Today, ".DotPlot.Targets.ColocGrouped.TCellonly.GenderFemales.ps"), plot = p2_females_cells)ggsave(paste0(PLOT_loc, "/", Today, ".DotPlot.Targets.ColocGrouped.TCellonly.GenderFemales.pdf"), plot = p2_females_cells)
# ggsave(paste0(PLOT_loc, "/", Today, ".DotPlot.Targets.ColocGrouped.TCellonly.Gender.png"), plot = combined_plot_cells)
# ggsave(paste0(PLOT_loc, "/", Today, ".DotPlot.Targets.ColocGrouped.TCellonly.Gender.ps"), plot = combined_plot_cells)
# ggsave(paste0(PLOT_loc, "/", Today, ".DotPlot.Targets.ColocGrouped.TCellonly.Gender.pdf"), plot = combined_plot_cells)FeaturePlot(scRNAseqDataCEA39, features = c(target_genes_qc),
cols = c("#ECECEC", "#DB003F", "#9A3480","#1290D9"),
combine = TRUE)
ggsave(paste0(PLOT_loc, "/", Today, ".FeaturePlot.Targets.png"), plot = last_plot())# VlnPlot(scRNAseqDataCEA39, features = "DUSP27")
# VlnPlot files
ifelse(!dir.exists(file.path(PLOT_loc, "/VlnPlot")),
dir.create(file.path(PLOT_loc, "/VlnPlot")),
FALSE)[1] FALSE
VlnPlot_loc = paste0(PLOT_loc, "/VlnPlot")
for (GENE in target_genes_qc){
print(paste0("Projecting the expression of ", GENE, "."))
vp1 <- VlnPlot(scRNAseqDataCEA39, features = GENE) +
xlab("cell communities") +
ylab(bquote("normalized expression")) +
theme(axis.title.x = element_text(color = "#000000", size = 14, face = "bold"),
axis.title.y = element_text(color = "#000000", size = 14, face = "bold"),
legend.position = "none")
ggsave(paste0(VlnPlot_loc, "/", Today, ".VlnPlot.",GENE,".png"), plot = last_plot())
ggsave(paste0(VlnPlot_loc, "/", Today, ".VlnPlot.",GENE,".ps"), plot = last_plot())
ggsave(paste0(VlnPlot_loc, "/", Today, ".VlnPlot.",GENE,".pdf"), plot = last_plot())
# print(vp1)
}[1] "Projecting the expression of SCGB3A2."
[1] "Projecting the expression of LIX1."
[1] "Projecting the expression of IGSF9B."
[1] "Projecting the expression of ALB."
[1] "Projecting the expression of OR10A4."
[1] "Projecting the expression of RASEF."
[1] "Projecting the expression of NEDD4."
[1] "Projecting the expression of TCL1A."
[1] "Projecting the expression of FBXO15."
[1] "Projecting the expression of F5."
[1] "Projecting the expression of TMEM212."
[1] "Projecting the expression of PTPRD."
[1] "Projecting the expression of CYP46A1."
[1] "Projecting the expression of OR2T33."
[1] "Projecting the expression of SORBS2."
[1] "Projecting the expression of ITGA7."
[1] "Projecting the expression of FOS."
[1] "Projecting the expression of HMOX1."
[1] "Projecting the expression of LAPTM5."
[1] "Projecting the expression of MMP9."
[1] "Projecting the expression of SIGLEC1."
[1] "Projecting the expression of FTL."
[1] "Projecting the expression of CD14."
[1] "Projecting the expression of HCST."
[1] "Projecting the expression of TIMP3."
[1] "Projecting the expression of CCL2."
[1] "Projecting the expression of SAT1."
[1] "Projecting the expression of CD163."
[1] "Projecting the expression of PTGDS."
[1] "Projecting the expression of LGALS9."
[1] "Projecting the expression of ACKR1."
[1] "Projecting the expression of NT5DC2."
[1] "Projecting the expression of TGFBI."
[1] "Projecting the expression of C1QC."
[1] "Projecting the expression of S100A9."
Here we project genes to only the broad cell communities:
[1] CD3+ TC I CD3+ TC IV CD34+ EC I
[4] CD3+ TC V CD3+CD56+ NK II CD3+ TC VI
[7] CD68+IL18+TLR4+TREM2+ MRes CD3+CD56+ NK I ACTA2+ SMC
[10] CD3+ TC II FOXP3+ TC CD34+ EC II
[13] CD3+ TC III CD68+CD1C+ DC CD68+CASP1+IL1B+SELL MInf
[16] CD79A+ BCmem CD68+ABCA1+OLR1+TREM2+ FC CD68+KIT+ MC
[19] CD68+CD4+ Mono CD79+ BCplasma
20 Levels: CD68+CD4+ Mono CD68+IL18+TLR4+TREM2+ MRes CD68+CD1C+ DC ... CD79+ BCplasma
Comparison between the macrophages cell communities (CD14/CD68+), and all other communities.
MAC.markers <- FindMarkers(object = scRNAseqDataCEA39,
ident.1 = c("CD68+CASP1+IL1B+SELL MInf",
"CD68+CD1C+ DC",
"CD68+CD4+ Mono",
"CD68+IL18+TLR4+TREM2+ MRes",
"CD68+ABCA1+OLR1+TREM2+ FC"),
ident.2 = c(#"CD68+CASP1+IL1B+SELL MInf",
#"CD68+CD1C+ DC",
#"CD68+CD4+ Mono",
#"CD68+IL18+TLR4+TREM2+ MRes",
#"CD68+ABCA1+OLR1+TREM2+ FC",
"CD3+ TC I",
"CD3+ TC II",
"CD3+ TC III",
"CD3+ TC IV",
"CD3+ TC V",
"CD3+ TC VI",
"FOXP3+ TC",
"CD34+ EC I",
"CD34+ EC II",
"ACTA2+ SMC",
"CD3+CD56+ NK I",
"CD3+CD56+ NK II",
"CD68+KIT+ MC",
"CD79+ BCplasma",
"CD79A+ BCmem"))
DT::datatable(MAC.markers)MAC_Volcano_TargetsA = EnhancedVolcano(MAC.markers,
lab = rownames(MAC.markers),
x = "avg_log2FC",
y = "p_val_adj",
selectLab = target_genes_qc,
axisLabSize = 12,
xlab = "average fold-change",
title = "Macrophage markers\n(Macrophage communities vs the rest)",
titleLabSize = 14,
pCutoff = 0.05/(nrow(MAC.markers)), # 20552 genes
FCcutoff = 1.25,
pointSize = 1.5,
labSize = 3.0,
legendLabels =c('NS','avg. fold-change','P',
'P & avg. fold-change'),
legendPosition = "right",
legendLabSize = 10,
legendIconSize = 3.0,
drawConnectors = TRUE,
widthConnectors = 0.2,
colConnectors = "#595A5C",
gridlines.major = FALSE,
gridlines.minor = FALSE)
MAC_Volcano_TargetsA
ggsave(paste0(PLOT_loc, "/", Today, ".Volcano.MAC.DEG.Targets.pdf"),
plot = MAC_Volcano_TargetsA)The target results are given below and written to a file.
Comparison between the smooth muscle cell communities (ACTA2+), and all other communities.
SMC.markers <- FindMarkers(object = scRNAseqDataCEA39,
ident.1 = c("ACTA2+ SMC"),
ident.2 = c("CD68+CASP1+IL1B+SELL MInf",
"CD68+CD1C+ DC",
"CD68+CD4+ Mono",
"CD68+IL18+TLR4+TREM2+ MRes",
"CD68+ABCA1+OLR1+TREM2+ FC",
"CD3+ TC I",
"CD3+ TC II",
"CD3+ TC III",
"CD3+ TC IV",
"CD3+ TC V",
"CD3+ TC VI",
"FOXP3+ TC",
"CD34+ EC I",
"CD34+ EC II",
#"ACTA2+ SMC",
"CD3+CD56+ NK I",
"CD3+CD56+ NK II",
"CD68+KIT+ MC",
"CD79+ BCplasma",
"CD79A+ BCmem"))
DT::datatable(SMC.markers)SMC_Volcano_TargetsA = EnhancedVolcano(SMC.markers,
lab = rownames(SMC.markers),
x = "avg_log2FC",
y = "p_val_adj",
selectLab = target_genes_qc,
axisLabSize = 12,
xlab = "average fold-change",
title = "SMC markers\n(SMC communities vs the rest)",
titleLabSize = 14,
pCutoff = 0.05/(nrow(SMC.markers)), # 20552 genes
FCcutoff = 1.25,
pointSize = 1.5,
labSize = 3.0,
legendLabels =c('NS','avg. fold-change','P',
'P & avg. fold-change'),
legendPosition = "right",
legendLabSize = 10,
legendIconSize = 3.0,
drawConnectors = TRUE,
widthConnectors = 0.2,
colConnectors = "#595A5C",
gridlines.major = FALSE,
gridlines.minor = FALSE)
SMC_Volcano_TargetsA
ggsave(paste0(PLOT_loc, "/", Today, ".Volcano.SMC.DEG.Targets.pdf"),
plot = SMC_Volcano_TargetsA)The target results are given below and written to a file.
Comparison between the endothelial cell communities (CD34+), and all other communities.
EC.markers <- FindMarkers(object = scRNAseqDataCEA39,
ident.1 = c("CD34+ EC I",
"CD34+ EC II"),
ident.2 = c("CD68+CASP1+IL1B+SELL MInf",
"CD68+CD1C+ DC",
"CD68+CD4+ Mono",
"CD68+IL18+TLR4+TREM2+ MRes",
"CD68+ABCA1+OLR1+TREM2+ FC",
"CD3+ TC I",
"CD3+ TC II",
"CD3+ TC III",
"CD3+ TC IV",
"CD3+ TC V",
"CD3+ TC VI",
"FOXP3+ TC",
# "CD34+ EC I",
# "CD34+ EC II",
"ACTA2+ SMC",
"CD3+CD56+ NK I",
"CD3+CD56+ NK II",
"CD68+KIT+ MC",
"CD79+ BCplasma",
"CD79A+ BCmem"))
DT::datatable(EC.markers)EC_Volcano_TargetsA = EnhancedVolcano(EC.markers,
lab = rownames(EC.markers),
x = "avg_log2FC",
y = "p_val_adj",
selectLab = target_genes_qc,
axisLabSize = 12,
xlab = "average fold-change",
title = "Endothelial cell markers\n(EC communities vs the rest)",
titleLabSize = 14,
pCutoff = 0.05/(nrow(EC.markers)), # 20552 genes
FCcutoff = 1.25,
pointSize = 1.5,
labSize = 3.0,
legendLabels =c('NS','avg. fold-change','P',
'P & avg. fold-change'),
legendPosition = "right",
legendLabSize = 10,
legendIconSize = 3.0,
drawConnectors = TRUE,
widthConnectors = 0.2,
colConnectors = "#595A5C",
gridlines.major = FALSE,
gridlines.minor = FALSE)
EC_Volcano_TargetsA
ggsave(paste0(PLOT_loc, "/", Today, ".Volcano.EC.DEG.Targets.pdf"),
plot = EC_Volcano_TargetsA)The target results are given below and written to a file.
Comparison between the T-cell communities (CD3/CD4/CD8+), and all other communities.
TC.markers <- FindMarkers(object = scRNAseqDataCEA39,
ident.1 = c("CD3+ TC I",
"CD3+ TC II",
"CD3+ TC III",
"CD3+ TC IV",
"CD3+ TC V",
"CD3+ TC VI",
"FOXP3+ TC"),
ident.2 = c("CD68+CASP1+IL1B+SELL MInf",
"CD68+CD1C+ DC",
"CD68+CD4+ Mono",
"CD68+IL18+TLR4+TREM2+ MRes",
"CD68+ABCA1+OLR1+TREM2+ FC",
# "CD3+ TC I",
# "CD3+ TC II",
# "CD3+ TC III",
# "CD3+ TC IV",
# "CD3+ TC V",
# "CD3+ TC VI",
# "FOXP3+ TC",
"CD34+ EC I",
"CD34+ EC II",
"ACTA2+ SMC",
"CD3+CD56+ NK I",
"CD3+CD56+ NK II",
"CD68+KIT+ MC",
"CD79+ BCplasma",
"CD79A+ BCmem"))
DT::datatable(TC.markers)TC_Volcano_TargetsA = EnhancedVolcano(TC.markers,
lab = rownames(TC.markers),
x = "avg_log2FC",
y = "p_val_adj",
selectLab = target_genes_qc,
axisLabSize = 12,
xlab = "average fold-change",
title = "T-cell markers\n(T-cell communities vs the rest)",
titleLabSize = 14,
pCutoff = 0.05/nrow(TC.markers), # 20552 genes
FCcutoff = 1.25,
pointSize = 1.5,
labSize = 3.0,
legendLabels =c('NS','avg. fold-change','P',
'P & avg. fold-change'),
legendPosition = "right",
legendLabSize = 10,
legendIconSize = 3.0,
drawConnectors = TRUE,
widthConnectors = 0.2,
colConnectors = "#595A5C",
gridlines.major = FALSE,
gridlines.minor = FALSE)
TC_Volcano_TargetsA
ggsave(paste0(PLOT_loc, "/", Today, ".Volcano.TC.DEG.Targets.pdf"),
plot = TC_Volcano_TargetsA)The target results are given below and written to a file.
Comparison between the B-cell communities (CD79A+), and all other communities.
BC.markers <- FindMarkers(object = scRNAseqDataCEA39,
ident.1 = c("CD79+ BCplasma",
"CD79A+ BCmem"),
ident.2 = c("CD68+CASP1+IL1B+SELL MInf",
"CD68+CD1C+ DC",
"CD68+CD4+ Mono",
"CD68+IL18+TLR4+TREM2+ MRes",
"CD68+ABCA1+OLR1+TREM2+ FC",
"CD3+ TC I",
"CD3+ TC II",
"CD3+ TC III",
"CD3+ TC IV",
"CD3+ TC V",
"CD3+ TC VI",
"FOXP3+ TC",
"CD34+ EC I",
"CD34+ EC II",
"ACTA2+ SMC",
"CD3+CD56+ NK I",
"CD3+CD56+ NK II",
"CD68+KIT+ MC"
# "CD79+ BCplasma",
# "CD79A+ BCmem"
))
DT::datatable(BC.markers)BC_Volcano_TargetsA = EnhancedVolcano(BC.markers,
lab = rownames(BC.markers),
x = "avg_log2FC",
y = "p_val_adj",
selectLab = target_genes_qc,
axisLabSize = 12,
xlab = "average fold-change",
title = "B-cell markers\n(B-cell communities vs the rest)",
titleLabSize = 14,
pCutoff = 0.05/nrow(BC.markers), # 20552 genes
FCcutoff = 1.25,
pointSize = 1.5,
labSize = 3.0,
legendLabels =c('NS','avg. fold-change','P',
'P & avg. fold-change'),
legendPosition = "right",
legendLabSize = 10,
legendIconSize = 3.0,
drawConnectors = TRUE,
widthConnectors = 0.2,
colConnectors = "#595A5C",
gridlines.major = FALSE,
gridlines.minor = FALSE)
BC_Volcano_TargetsA
ggsave(paste0(PLOT_loc, "/", Today, ".Volcano.BC.DEG.Targets.pdf"),
plot = BC_Volcano_TargetsA)The target results are given below and written to a file.
Comparison between the mast cell communities (KIT+), and all other communities.
MC.markers <- FindMarkers(object = scRNAseqDataCEA39,
ident.1 = c("CD68+KIT+ MC"),
ident.2 = c("CD68+CASP1+IL1B+SELL MInf",
"CD68+CD1C+ DC",
"CD68+CD4+ Mono",
"CD68+IL18+TLR4+TREM2+ MRes",
"CD68+ABCA1+OLR1+TREM2+ FC",
"CD3+ TC I",
"CD3+ TC II",
"CD3+ TC III",
"CD3+ TC IV",
"CD3+ TC V",
"CD3+ TC VI",
"FOXP3+ TC",
"CD34+ EC I",
"CD34+ EC II",
"ACTA2+ SMC",
"CD3+CD56+ NK I",
"CD3+CD56+ NK II",
# "CD68+KIT+ MC",
"CD79+ BCplasma",
"CD79A+ BCmem"))
DT::datatable(MC.markers)MC_Volcano_TargetsA = EnhancedVolcano(MC.markers,
lab = rownames(MC.markers),
x = "avg_log2FC",
y = "p_val_adj",
selectLab = target_genes_qc,
axisLabSize = 12,
xlab = "average fold-change",
title = "Mast cell markers\n(Mast cell communities vs the rest)",
titleLabSize = 14,
pCutoff = 0.05/nrow(MC.markers), # 20552 genes
FCcutoff = 1.25,
pointSize = 1.5,
labSize = 3.0,
legendLabels =c('NS','avg. fold-change','P',
'P & avg. fold-change'),
legendPosition = "right",
legendLabSize = 10,
legendIconSize = 3.0,
drawConnectors = TRUE,
widthConnectors = 0.2,
colConnectors = "#595A5C",
gridlines.major = FALSE,
gridlines.minor = FALSE)
MC_Volcano_TargetsA
ggsave(paste0(PLOT_loc, "/", Today, ".Volcano.MC.DEG.Targets.pdf"),
plot = MC_Volcano_TargetsA)The target results are given below and written to a file.
Comparison between the natural killer cell communities (NCAM1+), and all other communities.
NK.markers <- FindMarkers(object = scRNAseqDataCEA39,
ident.1 = c("CD3+CD56+ NK I",
"CD3+CD56+ NK II"),
ident.2 = c("CD68+CASP1+IL1B+SELL MInf",
"CD68+CD1C+ DC",
"CD68+CD4+ Mono",
"CD68+IL18+TLR4+TREM2+ MRes",
"CD68+ABCA1+OLR1+TREM2+ FC",
"CD3+ TC I",
"CD3+ TC II",
"CD3+ TC III",
"CD3+ TC IV",
"CD3+ TC V",
"CD3+ TC VI",
"FOXP3+ TC",
"CD34+ EC I",
"CD34+ EC II",
"ACTA2+ SMC",
# "CD3+CD56+ NK I",
# "CD3+CD56+ NK II",
"CD68+KIT+ MC",
"CD79+ BCplasma",
"CD79A+ BCmem"))
DT::datatable(NK.markers)NK_Volcano_TargetsA = EnhancedVolcano(NK.markers,
lab = rownames(NK.markers),
x = "avg_log2FC",
y = "p_val_adj",
selectLab = target_genes_qc,
axisLabSize = 12,
xlab = "average fold-change",
title = "NK markers\n(NK-cell communities vs the rest)",
titleLabSize = 14,
pCutoff = 0.05/nrow(NK.markers), # 20552 genes
FCcutoff = 1.25,
pointSize = 1.5,
labSize = 3.0,
legendLabels =c('NS','avg. fold-change','P',
'P & avg. fold-change'),
legendPosition = "right",
legendLabSize = 10,
legendIconSize = 3.0,
drawConnectors = TRUE,
widthConnectors = 0.2,
colConnectors = "#595A5C",
gridlines.major = FALSE,
gridlines.minor = FALSE)
NK_Volcano_TargetsA
ggsave(paste0(PLOT_loc, "/", Today, ".Volcano.NK.DEG.Targets.pdf"),
plot = NK_Volcano_TargetsA)The target results are given below and written to a file.
Let’s get the subset of most differentially expression genes for each subtype of macrophage, and compare these to the other cell communities.
library(tibble)
library(data.table) # for fwrite
# 1 - Define mapping of short prefixes to full cluster names:
cell_types <- c(
MAC_MINF = "CD68+CASP1+IL1B+SELL MInf",
MAC_DC = "CD68+CD1C+ DC",
MAC_MONO = "CD68+CD4+ Mono",
MAC_MRES = "CD68+IL18+TLR4+TREM2+ MRes",
MAC_FC = "CD68+ABCA1+OLR1+TREM2+ FC"
)
# 2 - Pull out the full list of clusters (so we can set ident.2 = everything else)
all_clusters <- unname(cell_types)
# 3 - Loop over each cell_type:
for (prefix in names(cell_types)) {
full_name <- cell_types[prefix]
# define ident.2 as “all clusters except the one in ident.1”
other_clusters <- setdiff(all_clusters, full_name)
# 4 - Run FindMarkers
markers <- FindMarkers(
object = scRNAseqDataCEA39,
ident.1 = full_name,
ident.2 = other_clusters
)
# 5 - add Gene column
markers <- add_column(markers, Gene = rownames(markers), .before = 1)
# 6 - Sort by p-value (most significant first)
markers_sorted_all <- markers[order(markers$p_val), , drop = FALSE]
top50_all <- head(markers_sorted_all, 50)
# 7 - Remove HLA- genes, then sort & take top 50
noHLA <- markers[!grepl("^HLA-", markers$Gene), , drop = FALSE]
markers_sorted_noHLA <- noHLA[order(noHLA$p_val), , drop = FALSE]
top50_noHLA <- head(markers_sorted_noHLA, 50)
# 8 - Write both files
base_name <- paste0(Today, ".", prefix)
fwrite(top50_all, file = file.path(OUT_loc, paste0(base_name, ".Top50All.txt")),
quote = FALSE, sep = "\t", showProgress = FALSE, verbose = FALSE)
fwrite(top50_noHLA, file = file.path(OUT_loc, paste0(base_name, ".Top50NoHLA.txt")),
quote = FALSE, sep = "\t", showProgress = FALSE, verbose = FALSE)
message("Wrote files for ", prefix, ":\n",
" • Top 50 (all genes): ", base_name, ".Top50All.txt\n",
" • Top 50 (no HLA-genes): ", base_name, ".Top50NoHLA.txt")
}Version: v1.2.3
Last update: 2025-07-02
Written by: Sander W. van der Laan (s.w.vanderlaan-2[at]umcutrecht.nl).
Description: Script to load single-cell RNA sequencing (scRNAseq) data, and perform quality control (QC), and initial mapping to cells.
Minimum requirements: R version 3.5.2 (2018-12-20) -- 'Eggshell Igloo', macOS Mojave (10.14.2).
**MoSCoW To-Do List**
The things we Must, Should, Could, and Would have given the time we have.
_M_
_S_
_C_
_W_
**Changes log**
* v1.2.3 Added a loop to extract the top 50 markers for each macrophage subtype.
* v1.2.2 Update to the writing of the baseline characteristics and tables.
* v1.2.1 Some fixes to the subsetting of a cell community and the visualization by group of target genes.
* v1.2.0 Fixed an issue where the subset of scRNAseq was loosing cell-identities.
* v1.1.1 Textual fixes.
* v1.1.1 Fix writing baseline table.
* v1.1.0 Update to study database.
* v1.0.2 Fixes to the start of the notebook. Update to loading of the clinical data. Fix on the gene-filtering.
* v1.0.1 Update to main AEDB (there is an error in the Age-variable in the new version). Fewer patients in scRNAseq (32 vs 39 with the newer dataset).
* v1.0.0 Initial version.
R version 4.5.1 (2025-06-13)
Platform: x86_64-apple-darwin20
Running under: macOS Sequoia 15.6
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: Europe/Amsterdam
tzcode source: internal
attached base packages:
[1] stats4 tools stats graphics grDevices utils datasets methods base
other attached packages:
[1] cowplot_1.1.3 RColorBrewer_1.1-3 flextable_0.9.7 officer_0.6.8
[5] writexl_1.5.4 gtsummary_2.2.0 labelled_2.14.1 openxlsx_4.2.8
[9] rmarkdown_2.29 worcs_0.1.18 Seurat_5.3.0 SeuratObject_5.1.0
[13] sp_2.2-0 devtools_2.4.5 usethis_3.1.0 BiocManager_1.30.25
[17] tableone_0.13.2 haven_2.5.4 EnhancedVolcano_1.26.0 ggrepel_0.9.6
[21] mygene_1.44.0 txdbmaker_1.4.1 GenomicFeatures_1.60.0 GenomicRanges_1.60.0
[25] GenomeInfoDb_1.44.0 org.Hs.eg.db_3.21.0 AnnotationDbi_1.70.0 IRanges_2.42.0
[29] S4Vectors_0.46.0 Biobase_2.68.0 BiocGenerics_0.54.0 generics_0.1.4
[33] DT_0.33 knitr_1.50 lubridate_1.9.4 forcats_1.0.0
[37] stringr_1.5.1 purrr_1.0.4 tibble_3.3.0 ggplot2_3.5.2
[41] tidyverse_2.0.0 data.table_1.17.2 naniar_1.1.0 tidylog_1.1.0
[45] tidyr_1.3.1 dplyr_1.1.4 optparse_1.7.5 readr_2.1.5
loaded via a namespace (and not attached):
[1] progress_1.2.3 urlchecker_1.0.1 nnet_7.3-20
[4] goftest_1.2-3 katex_1.5.0 Biostrings_2.76.0
[7] vctrs_0.6.5 spatstat.random_3.4-1 rticles_0.27
[10] digest_0.6.37 png_0.1-8 proxy_0.4-27
[13] deldir_2.0-4 parallelly_1.44.0 renv_1.1.4
[16] fontLiberation_0.1.0 MASS_7.3-65 reshape2_1.4.4
[19] httpuv_1.6.16 withr_3.0.2 ggrastr_1.0.2
[22] xfun_0.52 ellipsis_0.3.2 survival_3.8-3
[25] memoise_2.0.1 ggbeeswarm_0.7.2 profvis_0.4.0
[28] systemfonts_1.2.3 ragg_1.4.0 zoo_1.8-14
[31] V8_6.0.3 pbapply_1.7-2 Formula_1.2-5
[34] sys_3.4.3 prettyunits_1.2.0 KEGGREST_1.48.0
[37] promises_1.3.2 httr_1.4.7 restfulr_0.0.15
[40] globals_0.18.0 fitdistrplus_1.2-2 rstudioapi_0.17.1
[43] UCSC.utils_1.4.0 miniUI_0.1.2 base64enc_0.1-3
[46] curl_6.2.2 mitools_2.4 polyclip_1.10-7
[49] GenomeInfoDbData_1.2.14 SparseArray_1.8.0 xtable_1.8-4
[52] evaluate_1.0.3 S4Arrays_1.8.0 BiocFileCache_2.16.0
[55] hms_1.1.3 irlba_2.3.5.1 colorspace_2.1-1
[58] filelock_1.0.3 getopt_1.20.4 ROCR_1.0-11
[61] reticulate_1.42.0 spatstat.data_3.1-6 magrittr_2.0.3
[64] lmtest_0.9-40 glmGamPoi_1.20.0 later_1.4.2
[67] lattice_0.22-7 spatstat.geom_3.4-1 future.apply_1.11.3
[70] scattermore_1.2 XML_3.99-0.18 survey_4.4-2
[73] matrixStats_1.5.0 RcppAnnoy_0.0.22 class_7.3-23
[76] Hmisc_5.2-3 pillar_1.10.2 nlme_3.1-168
[79] beachmat_2.24.0 compiler_4.5.1 RSpectra_0.16-2
[82] stringi_1.8.7 tensor_1.5 SummarizedExperiment_1.38.1
[85] GenomicAlignments_1.44.0 plyr_1.8.9 crayon_1.5.3
[88] abind_1.4-8 BiocIO_1.18.0 chron_2.3-62
[91] bit_4.6.0 codetools_0.2-20 textshaping_1.0.1
[94] clisymbols_1.2.0 openssl_2.3.2 crosstalk_1.2.1
[97] bslib_0.9.0 e1071_1.7-16 plotly_4.10.4
[100] mime_0.13 splines_4.5.1 Rcpp_1.0.14
[103] fastDummies_1.7.5 sparseMatrixStats_1.20.0 dbplyr_2.5.0
[106] blob_1.2.4 fs_1.6.6 listenv_0.9.1
[109] checkmate_2.3.2 DelayedMatrixStats_1.30.0 prereg_0.6.0
[112] pkgbuild_1.4.7 gh_1.4.1 sqldf_0.4-11
[115] Matrix_1.7-3 statmod_1.5.0 tzdb_0.5.0
[118] visdat_0.6.0 equatags_0.2.1 pkgconfig_2.0.3
[121] cachem_1.1.0 RSQLite_2.3.11 viridisLite_0.4.2
[124] DBI_1.2.3 fastmap_1.2.0 scales_1.4.0
[127] grid_4.5.1 gt_1.0.0 credentials_2.0.2
[130] ica_1.0-3 cards_0.6.0 Rsamtools_2.24.0
[133] sass_0.4.10 broom_1.0.8 patchwork_1.3.1.9000
[136] cardx_0.2.4 dotCall64_1.2 RANN_2.6.2
[139] rpart_4.1.24 farver_2.1.2 gsubfn_0.7
[142] yaml_2.3.10 MatrixGenerics_1.20.0 foreign_0.8-90
[145] rtracklayer_1.68.0 cli_3.6.5 lifecycle_1.0.4
[148] askpass_1.2.1 uwot_0.2.3 presto_1.0.0
[151] sessioninfo_1.2.3 backports_1.5.0 BiocParallel_1.42.0
[154] timechange_0.3.0 gtable_0.3.6 rjson_0.2.23
[157] ggridges_0.5.6 progressr_0.15.1 limma_3.64.1
[160] parallel_4.5.1 jsonlite_2.0.0 RcppHNSW_0.6.0
[163] bitops_1.0-9 bit64_4.6.0-1 Rtsne_0.17
[166] spatstat.utils_3.1-4 proto_1.0.0 zip_2.3.3
[169] ranger_0.17.0 jquerylib_0.1.4 spatstat.univar_3.1-3
[172] lazyeval_0.2.2 shiny_1.10.0 xslt_1.5.1
[175] htmltools_0.5.8.1 sctransform_0.4.2 rappdirs_0.3.3
[178] tinytex_0.57 glue_1.8.0 spam_2.11-1
[181] httr2_1.1.2 XVector_0.48.0 gdtools_0.4.2
[184] RCurl_1.98-1.17 gridExtra_2.3 igraph_2.1.4
[187] R6_2.6.1 labeling_0.4.3 cluster_2.1.8.1
[190] pkgload_1.4.0 vipor_0.4.7 DelayedArray_0.34.1
[193] tidyselect_1.2.1 htmlTable_2.4.3 xml2_1.3.8
[196] fontBitstreamVera_0.1.1 future_1.49.0 KernSmooth_2.23-26
[199] fontquiver_0.2.1 htmlwidgets_1.6.4 biomaRt_2.64.0
[202] rlang_1.1.6 spatstat.sparse_3.1-0 gert_2.1.5
[205] spatstat.explore_3.4-2 uuid_1.2-1 remotes_2.5.0
[208] beeswarm_0.4.0
rm(backup.scRNAseqData)
rm(scRNAseqData, scRNAseqDataCEA39_females, scRNAseqDataCEA39_males)
# rm(scRNAseqDataCEA39) # Optional
save.image(paste0(PROJECT_loc, "/",Today,".",PROJECTNAME,".AESCRNA.results.RData"))| © 1979-2025 Sander W. van der Laan | s.w.vanderlaan[at]gmail[dot]com | vanderlaanand.science. |